/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Definitions of the full_sed_grid functions.

// $Id$

#include "full_sed_grid.h"
#include "grid-fits.h"
#include <CCfits/CCfits>
#include "blitz-fits.h"
#include <vector>
#include "density_generator.h"
#include "biniostream.h"
#include "fits-utilities.h"
#include "constants.h"
#include "grid.h"
#include "arepo_grid.h"

#ifdef WITH_AREPO
extern "C" {
#include "allvars.h"
#include "proto.h"
}
#endif

using namespace CCfits;
using namespace std;
using namespace blitz;

template<template<typename> class T>
mcrx::full_sed_grid<T>::full_sed_grid (boost::shared_ptr<T_grid_impl> g,
				       CCfits::ExtHDU& structure_hdu,
				       CCfits::ExtHDU& data_hdu,
				       const density_generator& gen,
				       int n_lambda) :
  polychromatic_grid<T> (g, data_hdu, gen, n_lambda)
{
  // this constructor only has to read units.

  structure_hdu.readKey("lengthunit", units_ ["length"]);
  Column& c_m_g = data_hdu.column("mass_gas");
  units_ ["mass"] = c_m_g.unit();
}


template<template<typename> class T>
mcrx::full_sed_grid<T>::full_sed_grid (boost::shared_ptr<T_grid_impl> g,
				       const T_unit_map& units):
  polychromatic_grid<T> (g)
{
  // this constructor only has to copy units.
  units_.insert(units.begin(), units.end());
  // however, this polychromatic_grid constructor doesn't call
  // make_blocks, so we need to do that here. The actual cells are
  // already set up in the grid, since we used a factory.
  this->make_blocks();
}


template<template<typename> class T>
mcrx::full_sed_grid<T>::full_sed_grid (boost::shared_ptr<T_grid_impl> g,
				       const density_generator& gen,
				       int n_lambda) :
  polychromatic_grid<T> (g)
{
#ifdef WITH_AREPO
  const int nc = this->n_cells ();
  assert(nc==N_gas);
  cout << "Creating cell data from Arepo data for " << nc << " cells. " << endl;

  // First copy units from the arepo units
  units_["length"]=arepo::arepo_units["length"];
  units_["time"]=arepo::arepo_units["time"];
  units_["mass"]=arepo::arepo_units["mass"];

  const T_float convert_density=
    All.UnitDensity_in_cgs*
    units::convert("g/cm^3",units_["mass"]+"/"+units_["length"]+"^3");
  const T_float convert_velocity=
    All.UnitVelocity_in_cm_per_s*
    units::convert("cm/s",units_["length"]+"/"+units_["time"]);

  int i = 0;
  bool negz=false;
  for (iterator c = this->begin (), e = this->end(); c != e; ++c, ++i) {
    // make an absorber with the appropriate densities. We only know
    // density of metals and gas so we ask the density generator to
    // give us a density vector, which we put into our big block.
    assert(c->volume()>0);
    const T_float rho_gas = SphP[i].Density*convert_density;

    // for some reason, Arepo gets tiny negative metallicities. If
    // this happens, set it to zero and print a warning.

    if(P[i].Metallicity<0) {
      negz=true;
      assert(P[i].Metallicity>-1e-6);
      P[i].Metallicity=0;
    }

    const T_float rho_met = rho_gas*P[i].Metallicity;
    vec3d vel=*reinterpret_cast<vec3d*>(&SphP[i].Momentum);
    vel/=(SphP[i].Density*SphP[i].Volume);

    const T_densities rho (gen.make_density_vector(rho_gas, rho_met));
    ASSERT_ALL(rho>=0);
    DEBUG(3,cout << "Cell " << i << " rho_gas=" << rho_gas << ", rho_met=" << rho_met << " rho=" << rho(0) << endl;);
    T_lambda intens (n_lambda);
    intens=0;
    T_data* d = new T_data (rho, vel, intens);
    c->set_data(d );
  }
  if(negz)
    cerr << "Warning: Negative metallicities found in Arepo grid.\n";

  // and now consolidate into continuous blocks
  this->make_blocks();

#else
  cerr << "Not compiled with Arepo support" << endl;
  exit(1);
#endif
}


/** Loads the radiation intensity in the grid
    cells from a FITS HDU. */
template<template<typename> class T>
void mcrx::full_sed_grid<T>::load_intensity (const CCfits::ExtHDU& i_hdu)
{
  cout << "Loading grid cell intensities from file" << endl;

  // Ideally, we would read directly into the big intensity array, but
  // because the values have to be converted using the cell volume, we
  // have to look anyway, so we just use a temporary array for
  // simplicity.

  // Indices are (scatterer, lambda).
  array_2 intensity;

  read(i_hdu, intensity);
  intensity = pow(10, intensity);

  float normalization = 1;
  T_float area_factor = 1;
  try {
    const_cast<ExtHDU&> (i_hdu).readKey ("normalization",
						 normalization);
  }
  catch (...) {}
  try {
    const_cast<ExtHDU&> (i_hdu).readKey ("area_factor",
						 area_factor);
  }
  catch (...) {}

  int ii=0;
  for (const_iterator c = this->begin (); c != this->end (); ++c, ++ii) {
    assert (ii < intensity.rows());
    // invert the calculation done in write_intensity here
    const T_lambda temp
      (intensity(ii, Range::all()) * 4*constants::pi*c->volume()*area_factor /
       normalization);
    // set intensity, remember this puts it into the big array.
    c->data()->get_absorber().set_intensity(temp);
  }
}

/** Dumps the radiation intensity to a file. */
template<template<typename> class T>
void mcrx::full_sed_grid<T>::write_dump (binofstream& file) const
{
  // We already have the intensities in one large array so this is easy.
  TinyVector<int, 2> extent = this->intensities_.shape();
  DEBUG(1,cout << "Dumping intensity data for grid, " << extent << " bytes" << endl;);
  file << extent [0] << extent [1];
  file.write(reinterpret_cast<const char*> (this->intensities_.dataFirst()),
	     this->intensities_.size()*sizeof (T_float) );
}

/** Loads the absorbed luminosity and radiation intensity in the grid
    cells from a previously written dump file. */
template<template<typename> class T>
void mcrx::full_sed_grid<T>::load_dump (binifstream& file) 
{
  // We should already have an appropriately sized large array for the
  // intensities, so loading is easy.
  TinyVector< int, 2> extent;
  file >> extent [0] >>  extent [1];
  DEBUG(1,cout << "Loading intensity data for grid, " << extent << " bytes" << endl;);
  assert(extent[0] == this->intensities_.shape()[0]);
  assert(extent[1] == this->intensities_.shape()[1]);

  file.read(reinterpret_cast<char*> (this->intensities_.dataFirst()),
	    this->intensities_.size()*sizeof (T_float) );
}


// Explicitly intantiate the templates necessary to avoid having to
// include everything

template class mcrx::full_sed_grid<mcrx::adaptive_grid>;
#ifdef WITH_AREPO
template class mcrx::full_sed_grid<mcrx::arepo_grid>;
#endif
